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AGN obscuration through dusty infrared dominated flows. II. 
Multidimensional, radiation-hydrodynamics modeling 
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ABSTRACT 

We explore a detailed model in which the active galactic nucleus (AGN) ob- 
scuration results from the extinction of AGN radiation in a global flow driven by 
the pressure of infrared radiation on dust grains. We assume that external illu- 
mination by UV and soft X-rays of the dusty gas located at approximately lpc 
away from the supermassive black hole is followed by a conversion of such radi- 
ation into IR. Using 2.5D, time-dependent radiation hydrodynamics simulations 
in a flux-limited diffusion approximation we find that the external illumination 
can support a geometrically thick obscuration via outflows driven by infrared ra- 
diation pressure in AGN with luminosities greater than 0.05 L e dd and Compton 
optical depth, r T > 1. 


1. Introduction 


A fundamental assumption of active galactic nuclei (AGN) unification schemes is that 
type 1 and type 2 AGNs have similar intrinsic properties. The basic premise of this paradigm 
is that obscuration and orientation effects are the major contributors to the observational 
dichotomy of AGNs. The goal of this paper is to suggest an approach which explains the 
AGN dichotomy as resulting from the extinction of the AGN radiation in a hydrodynamical 
outflow powered by the pressure of the infrared radiation on the dusty plasma of AGN 
outskirts. 


The suggestion that Seyfert 2 galaxies suffer from enhanced extinction compared to 
Seyfert 1 galaxies was made by Rowan-Robinson (1977) based on the infrared observations. 


However it was not until the seminal work of Antonucci (1984), and Antonucci & Miller 
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(1985) when key evidence was collected from studies based on optical spectropolarimetry. 
The detection of broad permitted lines in the polarized UV and optical spectrum of the 
nearby, luminous Seyfert 2 galaxy NGC 1068, confirmed that a bright, Seyfert 1 core is 
hidden behind optically thick, obscuring material. Notice that a prediction of polarization 
of the X-ray flux in the 0.1 — 10 keV range was made by Dorodnitsyn & Kallman (2010) 
based on theoretical modeling of AGN outflows. 

Direct evidence of the existence of toroidal obscuration comes from the mid- infrared 
observations of Seyfert 2 galaxies, such as the prototypical Seyfert 2 galaxy NGC 1068, 
and the closest AGN, the Circinus galaxy. Observations of NGC 1068 using VLTI reveal a 
multi-component, multi-temperature dusty conglomerate: an inner, relatively small ( 1 pc) 

and hot (~ 800 K) component embedded into an outer (~ 3.5 pc) component which is much 
colder (T ~ 320 K) ( Jaffe et al.|2004 Raban et al.||2009 ). In the Circinus galaxy the observed 
elongated, 0.4 pc in diameter component is interpreted as a disk-like structure seen almost 
edge-on. This disk-like structure is co-incident with that inferred from the VLBI maps of 


HoO maser emission (Greenhill et al. 2003), being embedded into a much larger rounded 


component. This is interpreted as a geometrically thick torus with temperature T < 300K 


(Tristram et al. 2007). 


The premise and the principal puzzle of AGN unification is the physical mechanism 
responsible for the geometrical thickness of the torus. Ample observational evidence for dust 
rules out support of the torus by gas pressure, as in such a case the temperature of the gas 
should be approximately of the order of the virial temperature, T v ir)g = 2.6 x 10 6 Mf/r pc K, 
where Mj is the black hole (BH) mass in 10 7 M 0 , and r pc is the distance in parsecs. 

Various mechanisms have been proposed to settle this issue: for example, in one of 


the first models a torus was 

considered being made of clumps having highly supersonic 

velocities ( 

Krolik & Begclman 

1988 

Beckert & Duschl 

2004 

). Magnetic fields are implicitly 


necessary in this model to provide enough elasticity to the clouds in order to avoid large 


dissipation though cloud-cloud collisions. Another model suggested by Phinney (1989), and 


Sanders et al. (1989) considered a locally geometrically thin, but globally warped disk. Global 


magnetic fields were suggested to be a key ingredient either though hydromagnetic winds 


(Konigl & Kartje 

1994 

), or as directly supporting the vertical balance of a quasi-static torus 

Lovelace et al. 

1998 

)• 

The main difficulty with hydromagnetic models comes from the large 


poloidal magnetic flux needed to support such a wind. It is unclear whether such a strong 
global poloidal magnetic field exists at large distances from a BH. The clumpy nature of 
an outflow was addressed by Elitzur & Shlosman (2006). These authors consider the dusty 


hydromagnetic obscuring wind as an alternative to quasi-static torus models. 


It has been pointed out by Pier & Krolik (1992) that the infrared radiation pressure 
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on dust may suffice to balance the vertical gravitational force. Based on these ideas (Krolik 


(2007), Shi & Krolik (2008) constructed a semi-analytic model of a static, infrared-supported 


torus. In all models in which infrared radiation is responsible for the torus thickness it is 
tacitly assumed that external radiation ranging from UV to soft X-rays is absorbed and 
converted into IR at the inner face of the torus. 

To sum up, previous models of AGN torus obscuration divide with respect to whether 
they are i) static, i.e. such as in a model of self-gravitating clouds or in a model of a static 
IR supported torus; ii) dynamic, such as in hydromagnetic wind scenario. 

In this paper we present simulations which demonstrate that obscuration at parses 
scale can be produced by global outflows driven by infrared radiation pressure on dust. 
Arguments for this model are both observational and theoretical. Simple estimates show 
that the observed dust temperatures (see above) translate into high infrared radiation energy 
densities. The latter when coupled with the high opacity of dust to IR radiation will produce 
strong radiation pressure force. Furthermore, one can easily see that radiation pressure on 
dust exceeds the gas pressure and that together with gravity and centrifugal forces they 


determine the fate of the torus. To support this line of arguments, in (Dorodnitsyn et al. 


(2011), Paper I), it was shown that if the temperature in the torus exceeds that of T v i r , r ~ 


312 (n/10 5 M 7 r pc 1 ) 1 / 4 — 987 (n/10 7 M 7 r pc 1 ) 1 / 4 K, where n is the number density, Mbh = 


pc 


1O 7 M 0 is the black hole mass, equilibrium between rotational, gravitational and radiation 
forces cannot be maintained resulting in the beginning of an outflow. Paper I also presented 
a numerical solution of 2D transfer and dynamics subject to the restriction that the outflow is 
ID and vertical. In Paper I it was assumed that external illumination by UV and soft X-rays 
and the subsequent conversion of this radiation into IR results in pumping of the torus with 
IR photons and producing a significant infrared pressure on dust. A conservative scenario 
for the radiation acceleration was employed, that is only external sources were considered 
with no contribution from the accretion disk. It was also assumed that a thin accretion disk 
provides the necessary mass loading for the wind. Such disk is buried in the dusty wind 
and its characteristic temperature can be of the order of a fewxlOO — 1000K. Thus the disk 
contribution to the infrared pressure on dust can be significant. Given the highly speculative 
nature of estimates for the viscous transport in self-gravitating disks we do not take the disk 
contribution into account in our conservative estimate for the torus radiation driving in this 
paper. 


Czerny & Hryniewicz (2011) suggested the dusty wind as a possible origin for the low- 


ionization part of the broad line region. The effective temperature of the disk where the low 
ionization part of a broad line flow is formed in this picture is T e g- ~ 1000K > T v ir,r- Thus, 
their broad line flow is the innermost part of the our obscuring flow. At these small radii the 
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uplifted dusty flow is exposed to external heating, the dust evaporates, and the radiation 
force quenches which result in a failed wind scenario. The boundary between the two regions 
is roughly set by the dust sublimation radius ~ 0.4(L/10 45 )pc, where L(ergs -1 ) is the total 
luminosity. 

The nature of the AGN obscuration problem calls for multidimensional simulations and 
reduce the predictive power of ID modeling (which is extremely successful, for example, in 
stellar evolution calculations). Contrary to Paper I, in the current work we solve the full 
system of time-dependent equations of radiation hydrodynamics in three dimensions with 
axisymmetry (2.5D), and find the fate of the dusty gas for conditions relevant to real AGNs. 

The plan of this paper is as follows: In Section [2] we discuss some of the basic physical 
properties of dust and gas in the torus, and review some of the results of Paper I which 
are extensively used in this paper. The torus problem can only be solved using radiation 
hydrodynamics (RHD). The basic setup of the equations of radiation hydrodynamics is 
explained in Section [3] We also describe how physical conditions in a radiationally dominated 
plasma of a torus influence various regimes at which such system of equations can be solved. 
In Section [4] we describe our numerical approach which we adopt in order to solve the 
full system of RHD. The implementation of the boundary conditions is also discussed. A 
description of the results obtained from a calculated grid of models is given in Section [5j We 
conclude with the discussion of the results and observational perspectives of our model in 
Section |6l 


2. Radiation and matter in the torus 

At the inner parts of the accretion disk near the BH copious UV-, and X-ray radiation 
is generated. The very high opacity of dusty plasma to UV radiation makes impossible 
any static spherically-symmetric configuration if the UV luminosity, L^ ust approaches the 
critical value (see Paper I): 


-^cldust — 5 x to -4 — 0.01 L edd , 
where L e dd is the Eddington critical luminosity: 


( 1 ) 


ToHd — 


A-kcGM} 


BH 


= 1.26 x 10 45 M 7 




( 2 ) 


where kt = 0.4 cm 2 g 1 is the Thomson opacity due to electron scattering, and M 7 = 
A/bh/( 10 7 M o ). Notice, that 0 is calculated taking k uv ~ 6 x 10 3 /vr, the dust grain sizes 
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of 0.025 — 0.25 //rn (Mathis et al.J[1977), grain density, n d = 2 — 3gcm , and 50 — 100 dust 
to gas mass ratio and assuming a perfect dust-gas coupling. For simplicity in the following 
k denotes the IR opacity of dust. 


Analyzing the simplified model of the radiationally and rotationally supported torus, 
in Paper 1 it was found the following approximate condition for the temperature, T for the 
beginning of an outflow: 

To>T vir , r (r 0 )ry 4 , (3) 

where Tq is the temperature at the base of the wind at the distance r o from BH, and 


/GW) 1/4 

\ ar J 

where p is the density, a = 7.56 • 10 -15 ergK _4 cm -3 is the radiation density constant, T c = 
L/L c , and L c is the critical luminosity in the infrared with respect to absorption on dust: 


4"7rc(jr 117 bH 

L c = — ^ 0.03 - 0.1 L edd , 

K 


(5) 


calculated assuming that the Rosseland mean opacity, k ~ 10 — 30kt in the temperature 
range 10 2 — 10 3 K (Semenov et al. 2003). 

The critical temperature Q is, in fact, a definition of the virial temperature in a 
radiation-dominated plasma. For densities relevant to numerical solutions presented in this 
paper T virjr spans from 312 (M 7 /r pc ) 1 / 4 K for number density, n = 10 5 , to 987 (M 7 /r pc ) 1 / 4 K 
for n — 10' . 


In the dusty plasma of a torus, at typical values of T and p the pressure is domi- 
nated by the radiation. Its relative importance is described by the parameter /3 = P g /P — 
(lO 3 T 3 /n-j + l) , where P is the total pressure 


p = p s + n, (6) 

where 

P g = — pUT, II = a T 4 /3, (7) 

p m 

where P g , and II are the gas and radiation pressures respectively, 77 = 8.31 • 10 7 erg KA 1 g _1 is 
the universal gas constant, and p m is the mean molecular weight. Thus, at typical densities 
n = 10 6 — 10 7 cm -3 and temperatures T = 10 2 — 10 3 K one finds P g II. 
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When the approximate condition ([3]) is met, a radiatively-driven wind develops. The 
terminal velocity, v°° then can be estimated assuming spherical symmetry and neglecting P g 
in favor of the radiation pressure, Fj&k/c, obtaining 


v°° = yj 2GM/r 0 (T c - 1) ~ 293 ( M 7 r 0 . 5 «; 10 ) 0 ' 5 kms" 1 . (8) 

The characteristic temperature of the conversion layer, T e g was found in Paper I: 


Tiff = 4 ctT 


,GMi 


0 

kt ar z j 


1/4 


~ 463 


-£705^^7x1/4 


K, 


pc 


( 9 ) 


where a ~ 0.5 is the fraction of the incident X-ray flux re-emitted in the IR into the torus, 
and T = T/Tedd In the calculations presented in this paper we adopt the parameter T, 
instead of T eS , calculating the latter from (|9]). 


3. The model setup and basic equations 


In the frame of reference co-moving with the fluid, the description of the interaction 
between radiation and matter is simplest, free from aberration and Doppler effects. For 
example, only in such a frame the adopted emissivity and absorptivity of matter are isotropic 
and only in such a frame they correspond to those tabulated from laboratory experiments. 
As such, it is important to distinguish between reference frames when casting equations for 
matter and radiation. The equations of radiation hydrodynamics, to first order in v/c, can 
be formulated in the following form (Mihalas & Mihalas 1984): 


D t p 

D t y 



p V • v = 0, 

Vp + g rad - V<f>, 

P 

-pV • V - 47TXP B + cxeE, 


— V • F — Vv : P + 47txpR — cXeE, 


( 10 ) 

(11) 

( 12 ) 

(13) 


where quantities related to matter: p, p, and e are the material mass density, gas pressure 
and gas energy density, v is the velocity; quantities related to radiation are the frequency- 
integrated moments: E is the radiation energy density, F is the radiation flux, and P is the 
radiation pressure tensor; y P , xe are the Planck mean and energy mean absorption opacities 
(in cm -1 ), ; c is the speed of light, B = <jT 4 /tt is the Planck function, and a = acj 4 is the 
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Stefan-Boltzmann constant, T is the gas temperature; other notation include the convective 
derivative, D t = + v • V, and Vv : P in (fl2|) denotes the contraction ( djVi)P l K Notice 

that all the dependent variables in ( 10 )-( 13 ) are evaluated in the co-moving frame. 


Frequency- independent moments E, F. which appear in the above set of RHD equations 
are obtained by calculating angular moments from the frequency-integrated specific intensity, 
J(r, fi, v, t ): 


E( r,t) = - j dv j) dfl J(r,fi, v,t), 

/ oo r 

dv J) dQ n/(r, fi, v, t). 


(14) 

(15) 


The frequency- independent radiation pressure tensor P is found from 

P(r ,t) — — J dv <j> dQ nn /(r, Q, v, t). 

The radiation force, g ra( j is calculated from the following relation 


(16) 


§rad 


1 xf F 
P 


(17) 


where Xf = Xa + Xt is the total flux mean opacity consisting of absorption opacity, y a and 
the Thomson scattering opacity, xt ■ I n the following we will not differentiate between x.f, 
Xf and Xe, and omit subscript from x where appropriate. 


In addition to ( 10 )-( 13 ), the full system of equations of radiation hydrodynamics should 


include an equation for F . However, we adopt a flux-limited diffusion approximation (FLD), 
and there is no need of such equation as the closure relation between F and E is found from 
the diffusion law: 


where 


F = —VE = —D'VE, 

3 Kp 


(18) 


D = c A, (19) 

is the diffusion coefficient and A is the photon mean free path: A = 1 /(up), where k — x/p- 

The diffusion approximation is adopted by tacitly assuming that optical depth r > 1. 
To take into account the possibility of r < 1 the diffusion approximation should be modified 


so it has a correct limiting behavior . Notice that in a free-streaming limit |F| — > cE. 
However, when r 1 the mean free path, A — > oo, and D — > oo, and |F| — > oo. That is, 
when optical depth becomes small, or when p — > 0, the standard diffusion approximation is 
no longer applicable. 


In order to overcome this problem the standard approach is to adopt the flux-limited dif- 
fusion approximation (Alme & Wilson 1974 Minerbo 1978; Levermore & Pomraning 1981). 
In an FLD approximation A is replaced by A* = A A, where A is the flux limiter. The flux 
limiter we adopt in the current work is that of Levermore & Pomraning (1981): 


2 + Alp 
6 + 37?lp + 


( 20 ) 


where Alp = A \S7E\/E. If r — )■ 0, then Alp — » oo, and |F| ~ cE. In the optically thick 
limit Alp — > 0 and A — > 1/3. 


The gas is assumed to be in a local thermodynamic equilibrium (LTE) at a temperature, 
which can be different from that of the effective temperature of radiation. The opacities 
are treated in a grey approximation, i.e. frequency-independent, constant dust opacity is 
assumed. The equation of state is taken to be that of an ideal polytropic gas, p = (7 — l)e 
with the ratio of specific heats 7. The gas temperature is obtained from the relation T = 
(7 — l)p m e/(plZ). Given the great variety of physical conditions in dusty molecular gas, 
exposed to UV, and X-ray radiation, we set p m = 1 throughout this paper. 


Notice that in Paper I, radiation properties were calculated assuming stationarity and 
radiative equilibrium. To understand why V • F = 0 may be not too bad an approximation 
it is instructive to compare time scales of the variation of the radiation held. In most cases 
we expect the radiation held to follow matter with the huid-how time scale, tf = l/v ~ 
l/i\ ~ 5 • 10 3 If Mj 1 / 2 yr, where is the orbital velocity, l is a typical length scale of the 
system, and l\ = 1/ lpc. Another important scale is the time which is necessary for a photon 
to travel distance l : t T — l/c ~ 3.43 Zi yr, for an optically thin case. In an optically thick 
case a photon travels across l by means of a random walk. The corresponding time scale in 
a diffusion regime is the diffusion time t d = l 2 / c\ = t x lj A ~ 530 yr 

If A <C tf the radiation held adjusts almost instantaneously to changes of the physical 
conditions in the how. Consequently, calculations are simplified immensely because the ex- 
plicit time variation of the radiation held can be ignored. At a given time step the properties 
of the how can be calculated taking into account the radiation as it is frozen i.e. at the how 
time scale the radiation is essentially a sequence of snapshots which instantly adapts to the 
how. 
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If additionally, optical depth r is sufficiently high then the radiation-matter energy 
exchange term in (12), and (13) scales as B — cE/Att ~ 0(B/t 2 ). If this condition is 


augmented by a requirement of a strict stationarity: d/dt = 0, the resultant situation is 
described as a radiative equilibrium and the radiation field can be found from the equation 
V • F = 0. As was mentioned, the latter was used in Paper I to find radiation field in a 
stationary outflowing torus. 

Contrary to Paper I, in the present work we do not assume the existence of the wind 
nor are we concerned with a wind solution which is stationary in a strict sense (both are the 
assumptions of Paper I). Rather, we solve for the full time dependence of both the radiation 
field and the flow. The effective coupling outlined above between time scales and local optical 
depth (see e.g. Bisnovatyi-Kogan & Blinnikov (1978)) presents a significant challenge to any 


RHD simulations. Careful analysis of the RHD system of equations (Mihalas & Mihalas 1984) 


prescribes that in order to be consistent in all regimes of the radiation-matter interaction, 
including the ability to describe correctly regions of the flow where r < 1 or r > 1, all terms 


in the system of RHD (10)- (13) must be retained. 


4. Solution: method 


4.1. Nondimensionalization 


In order to explicitly extract governing parameters it is convenient to convert ( 10 )-( 13 ) 


into a dimensionless form adopting dimensionless variables: f = r/R 0 , t = t/to, where Ro , 
is a fiducial distance from the BH, to = Rq/vq is the characteristic flow-time, and v = v/vq, 
where v 0 = 0c/ 2 , and 0 O = GM B h/ Ro', matter variables convert as p — p/p 0 , p = p/e 0 , 
e = e/e 0 , where p 0 = n 0 m p is the fiducial mass density and n 0 is the number density, m p is 
the proton mass, and eo = PoVq] radiation variables transform as E = E/E 0 . F = F/cE 0 , and 
P = P/Eq, where Eo = aUg, where To is the fiducial temperature. The opacity transforms 
as X — x/xo, where Xo = l/7?o- Using such nondimensionalization and further simplifying 
notation by (hereafter) omitting tilde, we obtain 


D t p = -pV-v, (21) 

pD t v = — Vp + AixF — pV$, (22) 

PA (!) = -/^V.F-Vv : P + l K p(T 4 - S ), (23) 

PAQ = -pV-v-^A lK p(T 4 -F), (24) 
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where /3 0 = Vq /c, and non-dimensional parameters are 


aT 0 4 


A , = * = 

^■1 o ? 

<?o PoVo 


M = (7 - 1) 


^OPm 

PqR,Tq 


( 25 ) 


We solve equations (21 )-(24 ) adopting a time-dependent, axisymmetric 2.5 D approxi- 


mation, meaning that we keep track of all cp- components of all vector quantities, such as 
rotational velocity, v^, but assume d/d<p = 0. A cylindrical z, R coordinate system is adopted 
in all our computations. 


We derive our code from the family of the ZEUS codes (Stone & Norman 1992a). 


The hydrodynamical part of our code partially adopts methods and infrastructure found 
in ZEUS-2D, and ZEUS-MP codes (Hayes et al. 2006), including our original modifications 


summarized in Dorodnitsyn et al. (2008). Our original implementation of the radiation 


module which is adopted in the current work is built upon a simplified version developed in 
Paper I. 

Notice that the hydrodynamical part of the ZEUS codes is based on a two-step update 
of the dependent variables: In the source step the local update is done for <9 t v, to account 
for various forces: gas pressure, (Vp)/p, radiation pressure g rac j, and the gravitational force 
V</>; d t e is updated taking into account p ■ Vv tern, i.e. pdV work, where V is the specific 


volume; shocks are treated adopting the artificial viscosity prescription of von Neumann & 


Richtmyer (1950). Corresponding viscous stresses and dissipation due to artificial viscosity, 
are also added in the source step. Next is the transport step: the previously updated 
quantities are further updated taking into account fluid advection. 

According to the general strategy adopted in the ZEUS codes the update of the radiation 


energy is also made adopting the operator splitting of equation ( 23 ) into source and transport 


terms. The most time consuming part is the solution of the equation (23) for E in the source 


step: Thus, the radiation source step includes the finite difference update of E from 


BE 

— = -/^V • F = /V X V ■ (DVE) 


= VP 




(26) 


Equation (26) is solved numerically adopting an alternative direction implicit scheme (ADI) 
(Fletcherfl988; Fedoren ko|[r99 4) . The fact that the effective time step in (26) is ~ c/v times 


larger than the original time step demands the sub-cycling when advancing (26) over the 
original time-step dt. This is done by splitting dt, into 10-1000 sub-steps. 
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Since in most (but not all!) of the dense parts of the flow t T ,t d <C £ f one obtains that 
solving for the V • F ~ 0 term is the most important part of the solution for the radiation 
field at a given time step. A three-point finite differencing stencil is used to approximate 


the diffusion operator in (26), which involves a solution of a tri-diagonal matrix equation. 


Further details of how the V • F term in equation (26) is taken into account can be found in 
Paper I. 

Our radiation module adopts and enhances the radiation module described in Paper I. 
The latter was designed to integrate equation 0 = dE/dw = — V ■ F ~ 0 over a pseudo-time 
w towards a stationary solution. This extensively tested module was augmented by the 


algorithms needed to take into account Vv : P + (T 4 _ E) terms in equation (26 ). This 


( 2001 ) 


part of the solution implements the well tested algorithms found in Turner & Stone 
version of the ZEUS code. 

The numerical grid used in the radiation module adopts the same staggered grid used 
in the hydrodynamic part of the ZEUS code. On such a grid, for example, E is placed cell- 


centered, while D is face centered. Finite differencing in (26) is done by approximating the 
derivatives making use of volume elements dl 2 = xdx and dl\ = dz. This approach allows 


to avoid approximation errors near the coordinate singularities (Stone & Norman 1992b). 


An important part of the method is the implementation of the implicit time-stepping 
algorithm, and from this, the most important is that in our method during the update from 
E™j to F” +1 the diffusion coefficients Dij are taken at the ’’old time”, t n . The method we 


are using for the radiation-matter interaction part adopts that of Turner & Stone (2001). 


This includes a simultaneous solution of the following finite difference algebraic system of 


equations obtained from (23)- (24): 


E n+ 1 - E n = (-(Vv : P ) n+1 + /3i fiUp n ((T n+1 ) 4 - E n+1 )) St , (27) 

e n+ i — e n = (— (pV-v) n - A 1 p 1 K n p n ({T n+1 ) 4 -E n+1 ))6t, (28) 

(29) 


where St = t n+l —t n is the time step, and the omitted subscript ij is assumed for all dependent 
variables, and other notation is introduced: /3i = Pq 1 , T n+1 = (T” +1 ) = A 2 (e n+1 )/p n . 


Benefits from such a discretization of variables in time are twofold: i) no linearization of 
difference equations is required, as opposed to a fully implicit method; ii) it allows for a much 
simpler update of the radiation and gas energies at the radiation-matter interaction step, in 


which case equations (27)-(28) are reduced to a single scalar algebraic equation (Turner & 


Stone 2001). The numerical solution of the latter is much more robust compared to what is 
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needed in a fully implicit method (e.g. such as in Hayes et ah (2006)). 


In our hydrodynamic framework we combine the original radiation diffusion and bound- 
ary conditions module (aka the ’’radiation module”), the original matter- and radiation- 
energy update module and the hydrodynamical module. The hydrodynamical part is taken 
from that of ZEUS-MP framework without alteration. This includes radiation advection step 
(transport step in zeus terminology) and the time-step calculation machinery. This part of 


the ZEUS-MP has been extensively tested by Hayes et al. (2006) against the Marshak wave 


test and the radiation-dominated shock wave test. 

The radiation module developed in Paper I was verified against the time evolution of 


the solution with only the diffusion operator in (23), the interaction between radiation and 


matter via matter-coupling terms in (23), (24), and the verification of the interface between 


the radiation terms in the rest of the zeus hydrodynamical algorithms. 


The tests we have implemented are those described in detail in Turner & Stone (2001), 


and we outline them in the following. In the first test we follow how the static uniform 
matter initially out of balance with radiation approaches thermal equilibrium. Assuming 


that E is constant (notice that in the torus E e), equation (24) is reduced to an ordinary 
differential equation (ODE) for unknown e. Solving this ODE we compare the resultant 
solution with the full solution from the radiation module and thus verify the radiation- 
matter energy update in the code. We found an excellent agreement (|e co de — e a | /E & < 1%), 
where e a is the solution of the ODE. 

The second major test involves the evolution of the radiation flux divergence term alone. 
Generally, the easy way to do it is to compare the numerical solution in the optically thick 
regime with the analytical solution of the heat diffusion equation. The idea is that one can 
write an analytic solution (such as (Turner & Stone 2001), eq. (47)) of the diffusion equation 
with constant diffusion coefficients on the unit square with periodic boundary conditions, 
and to compare it with the solution from the diffusion solver (with no hydrodynamics). This 
was done in Paper I where a very good agreement was obtained between such an analytic 
solution and the numerical solution of the 2D diffusion equation. 


4.2. Boundary conditions and the grid 

The boundary conditions (BC) are adopted from Paper I where one can End a cor- 
responding discussion. Cylindrical {z, R} coordinates were adopted. The computational 
domain spans from R 0 to Ri in horizontal, and from Zq = 0 to Z\ in vertical direction in 
the meridional plane. The stiffness imposed by radiation terms together with the necessity 
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to sub-cycle during the radiation time-step, matter-, and radiation-update steps makes the 
computation numerically entensive. For a given set of input parameters, we solve our prob- 
lem on a 200 x 200 grid, and compare the results with those obtained on low resolution, 
100 x 100 grid. On both resolutions the solution converges towards a quasi-stationary one. 

The BC for the matter provide p, e and v as follows: At the equatorial plane the inflow 
BC are adopted, having always v Ztin > 0. An outflow, i.e. v z < 0 through the equatorial 
plane is not permitted. Notice that we allow for v Z) - m to be arbitrary small, being numerically 
limited to a fraction of v 8ig , the sonic velocity calculated using only P g . During numerical 
calculations the value of u 2) j n mirrors the value of the velocity found in the adjacent cell inside 
the computational domain. Thus, u 2iin adjusts in the process of computation. The difference 
of our approach from others often adopted in simulations of accretion disk winds is that we 
obtain v Z) \ n and thus the mass-loss rate, M self-consistently. A power-law distribution for 
the density is assumed at the equator: 


p(0, x) = x d , (30) 

At all other boundaries including the inner, left boundary, outflow BC are adopted. The 
azimuthal component of the velocity v $ is assumed to be Keplerian at the equatorial plane. 

For the radiation at the left boundary which is located at the distance P 0 from the BH, 
we specify the distribution of energy density 


E(z,x 0 ) = E x0 z e , (31) 

where E xo = E{z = 0, Rq)/E 0 . At the equatorial plane, the disk ” photospheric” conditions 
are mimicked using the effective temperature, T e g, which provides 


F z (zq,R) = —DdE/dz = aT e 4 f 


(32) 


where T eff is calculated from a self-consistent ’’photospheric” boundary condition, i.e. T eS = 
T(zo,R). The free-streaming boundary conditions are assumed at the upper boundary, 
located at z\ and at the right boundary at R\. 


Thus, equations ( 21 )-( 24 ) with boundary conditions (30)- (32) are integrated in time 
until a quasi-static solution is found. 
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5. Results 

One of the primary objectives of this paper is to prove by multidimensional radiation 
hydrodynamics calculations the suggestion made in Paper I, that the structure which is 
usually referred to as an AGN torus is better described in terms of a radiationally supported 
flow rather than as a quasi-static obscuring torus. 

To approach this goal a set of models is calculated from the first principles of radiation 
hydrodynamics. The structure we model is represented by an extended, rotating ring of 
radiatively-dominated plasma in the {z, R} plane in cylindrical coordinates. Our calculations 
are performed in 2.5D meaning that we are able to calculate arbitrary 2D, time-dependent 
distributions of gas and radiation and to model an arbitrary axially-symmetric (i.e. d $ = 0) 
velocity field. Apart from the fundamental properties of the torus such as the distribution 
of mass, radiation energy density and velocity, the mass-loss rate and the torus mass are 
readily obtained from such calculations. 


5.0.1. Basic parameters 


The BH mass, Mbh = 1 x 10 7 M 0 , and the size of the obscuration region Rq = lpc are 
fixed at these values in all our models. Also the models are characterized by the Thomson 
optical depth of the torus, which is calculated at the inclination 6 = 7t/2 from the z-axis: 
Tt = / 0 °° KxpdR, and by parameter T which measures the intensity of the illumination of the 
inner face of the torus by UV and soft X-rays. The effective temperature of the conversion 
layer, T e g is found then from ([9]), and we also adopt T 0 = T e s in (25). The parameter cl = 0.5 
enters the equatorial distribution of density (30), and parameters E x0 = 1, and e = 0.1 shape 
the distribution of E in (31) and are fixed in all models. 


The dust opacity k is fixed at a constant value throughout the computations. The 
latter approximation may create an artificial situation when a very low density wind with 
n <C no is accelerated to very high velocities. This wind is an artifact of the adopted 
approximation: it has a negligible column density and a negligible mass-flux, and in a more 
detailed calculation, when n is allowed to depend on density such wind would not exist. 
It is convenient to introduce the following definitions: the average bulk velocity of the flow 

(v) = / pvdV / / pdV , where V is the total volume occupied by the flow; and the maximum 

Jv Jv 

velocity of the dense wind: found throughout the flow given the condition p(z, R) > p t h 

is satisfied. Experimenting with various threshold values p t h we found that if p t h is in the 
range p t h = 0.01 — 0.001p 0 the value of i d aY remains almost unchanged. Notice that such 
a definition of u* nax gives results which are approximately in accord with estimates based 
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on the kinetic output of the wind u k ; n — 2 L kin /M, where L kin = J pv 3 /2dT, is the kinetic 
luminosity of the wind, and £ is the outer boundary of the computational domain. 


Model 

r 

Rq 

r T 


n 0 

{v) 

v * 

max 

A k in 

L 

bol 

M 

1 

0.5 

1 

0.53 

3 

• 10 5 

138.4 

666.61 

1.23 • 10 41 

6.24 

• 10 44 

1.71 

2 

0.3 

1 

0.53 

3 

• 10 5 

104.87 

560.27 

6.38 • 10 4 ° 

3.74 

• 10 44 

1.38 

3 

0.1 

1 

0.53 

3 

• 10 5 

55.82 

373.46 

1.48 • 10 4 ° 

1.24 

• 10 44 

1.85 

4 

0.05 

1 

0.53 

3 

• 10 5 

36.63 

282.23 

6 • 10 39 

6.25 

o 

CO 

0.64 

5 

0.8 

1 

1.8 

1 

• 10 6 

112.62 

765 

4.11 • 10 41 

9.99 

• 10 44 

5 

6 

0.5 

1 

1.8 

1 

• 10 6 

90.55 

513.42 

2.56 • 10 41 

6.24 

• 10 44 

4.18 

7 

0.3 

1 

1.8 

1 

• 10 6 

73.18 

348.45 

1.37- 10 41 

3.74 

• 10 44 

3.37 

8 

0.1 

1 

1.8 

1 

• 10 6 

59.73 

159.74 

3.09 • 10 4 ° 

1.24 

• 10 44 

1.88 

9 

0.05 

1 

1.8 

1 

• 10 6 

42.8 

129.39 

1.31 • 10 4 ° 

6.25 

o 

CO 

1.39 


Table 1. Models characterized by initial parameters: T, i? 0 (P c )> Tr, characteristic den- 
sity no(cm~ 3 ), and the resulting averaged flow velocity, (u)(kms _1 ), the averaged maximum 
velocity, u* iax (km s~ 1 ) , the kinetic and bolometric luminosities, L kin (ergs _1 ), Lbo^ergs^ 1 ), 
and mass-loss rates M(M 0 yr _1 ). 

Table 1 describes the set of calculated models, summarizes the governing parameters 
and outlines the most important results. The models are divided into two broad categories 
with respect to the total Thomson optical depth at the equator: marginally optically thin 
Tt — 0.53 and optically thick tt — 0.8. For simplicity, we refer to the first category as 
type I models and to the second as type II. After t ~ 3 to dynamical times, where f 0 — 
1.5 x 10 11 rpc 2 M^ 1 ^ 2 s, the quasi-equilibrium solutions are found for all models. Characteristic 
distributions of p, E and v for type I models are shown in Figure [l](3j and for type II models 
in Figure |4}[7j 


5.0.2. The kinetic energy of the wind 

The first apparent trend with respect to the energy budget is that pumping radiation 
energy into the flow (larger T) results in a larger kinetic output, L kin in such models. For 
example we have F k m/Tboi = 9.6 • 10 5 for Model 4, with T = 0.05 while for Model 1 with 
T = 0.5 we get L kin /Lboi = 2 • 10 -4 . 

Similar comparison for marginally optically thin models with similar T gives L kin / Lboi = 
2 • 10~ 4 for Model 9, and L kin /Lb 0 i = 4.1 • 10 -4 for Model 6. A note of caution: the mass-loss 
rate, M in Table 1 reflects only the formal integral over the mass-flux at the boundary of 
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the computational domain, which may not reflect the true escape of the gas to infinity. To 
illustrate this, in the following we show that most of the gas remains gravitationally bound 
to the BH. 


log(p) 
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Fig. 1. — Model 1. Color-intensity plots of the dimensionless density, p (left) and dimension- 
less infrared radiation energy density, E (right), and the superimposed velocity field. Axes: 
distance in parsecs. 

Figure 1 shows the color-intensity plot of p (with superimposed velocity field ) and E 
from Model 1. The effective temperature of the radiation decreases from its maximum at 
the left boundary relatively smoothly, while the distribution of p demonstrates a pronounced 
disk-like structure. 

This is one of the most important results of the current work: the self-consistent distri- 
bution of density naturally evolves into a geometrically thick disk-like structure. Notice that 
only the density BC at the equator and the radiation energy at one boundary are provided. 
The aspect ratio of the resultant toroidal structure is h/R ~ 1, where h(R) is the vertical 
extent of the disk. Formation of a geometrically thick disk is observed in other models which 
have higher tt: in Figure [4] for Model 6, and in Figure [5] for Model 7. 


17 



Fig. 2. — Model 3: the surface plot of the z— and R— velocity components, where U esc is 
the escape velocity. Horizontal: R: distance from the BH in parsecs; z: distance from the 
equatorial plane in parsecs; 



Fig. 3. — Model 3: Color-intensity plot and contours of the total velocity v = ( u 2 z + n|.)) 1//2 , 
where U esc is the escape velocity. Horizontal: R: distance from the BH in parsecs; z: distance 
from the equatorial plane in parsecs; 
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Fig. 4. — Model 6. Color-intensity plots of the dimensionless density, p (left) and dimension- 
less infrared radiation energy density, E (right), and the superimposed velocity held. Axes: 
distance in parsecs. 
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Fig. 5. — Model 7. Color-intensity plots of the dimensionless density, p (left) and dimension- 
less infrared radiation energy density, E (right), and the superimposed velocity held. Axes: 
distance in parsecs. 
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Fig. 6. — Model 9: the surface plot of the z— and R— velocity components, where U esc is 
the escape velocity. Horizontal: R : distance from the BH in parsecs; z: distance from the 
equatorial plane in parsecs; 



Fig. 7. — Model 9: Color-intensity plot and contours of the total velocity v, where U esc is 
the escape velocity. Horizontal: R: distance from the BH in parsecs; z: distance from the 
equatorial plane in parsecs; 
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5.0.3. The structure of the velocity field in the torus 


The analysis of the velocity field shown by arrows in Figures 1,3, and 4 supports the basic 
hypothesis made at the beginning of this paper: the torus is formed by gas which is not in a 
static equilibrium. It is not surprising that the characteristic density, Uq (or alternatively 7 t) 
is one of the most important driving parameters, which critically influences the distribution 
of the radiation energy density and to a large extent determines the topology and magnitude 
of the velocity. The most important observation is that the disk-like structure seen in density 
plots represents a portion of a wind which is mostly subsonic. This is a striking analogy 
with models of outflowing stellar atmospheres driven by radiation pressure in the continuum 
(Bisnovatyi- Kogan & Dorodnit syn|[l999 ) where the inner subsonic part of the wind is usually 
the most dense one and the transition to r < 1 happens soon after the sonic point is reached. 


Figures [3j and [7] show the color- intensity plot and contours of the total velocity v = 
(u 2 z + ujf)) 11 ' 2 . One can see that in the first case the geometrically thick obscuring flow does 
not have enough speed to escape the potential well of the BH. The fast component is formed 
only when density drops. Low luminosity model shown in Figure [7] has a weak outflow 
without such a clear separation into low velocity and dense wind (i.e. ’’atmosphere” in a 
stellar analogy) and the fast wind. 

The difference between the lower density models (tt = 0.53, i.e. Figure [I]) and higher 
density ones (tt = 1.8, i.e. Figure [5]) is that in the latter there is a larger portion of the wind 
where the velocity is comparable to the escape velocity of the gas, U esc = (GM B h/ Rq) 1 ^ 2 — 
207 My^ 2 i7pc^ 2 kms _1 . Most of the torus mass is participating in the intensive global motion, 
though most of the gas remains in the potential well of the BH. Notice that we do not see 
a quasi-static disk in our solution at any stage of our simulation. However if T is not too 
small the velocities in a denser, disk-like part of the flow are such that vr <C v z ~ U esc (c.f. 
Figure [2]). A similar situation is observed in the most of our models except for type II models 
with T = 0.05. In the latter case (c.f. Figure [6]), vr ~ v z everywhere in the computational 
domain, except for the region close to the conversion layer, where v z > vr. The results show 
that the average velocity of the obscuring flow, ( v ) is almost always comparable but smaller 
than the escape velocity: for example, (v) = 0.53 U esc for Model 5; ( v ) = 0.43 U esc for the 
Model 6; and (v) = 0.28U esc for Model 8. 


For models with tt = 0.53 the maximum velocity reaches v* iax = 3.17 U esc for Model 
1) ?, max = 2.66 7/ eS c for Model 2, dropping further until it is n^ ax = 1.34 U esc for Model 4. 
Comparable results are obtained for models with tt = 1.8: the maximum nj^ ax = 3.64 U esc 
is obtained for Model 5. A T is reduced so does the being n^ ax = 2.44 7/ esc for Model 
6, and finally r?* iax = 0.76 U esc is obtained for Model 8. 
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5.0.4 ■ The column density and optical depth 

Comparing various models at 0 — 45° one can see that Models 1-4 (tt = 0.53) do not 
have enough column density to provide Compton-thick obscuration at this inclination. On 
the other hand, at higher inclinations both sets of models provide a noticeable extinction of 
the light from the supermassive BH. Notice that we do not account for the possible nonzero 
optical depth provided at smaller radii by any sort of additional wind, such as MHD wind 
or warm absorber flow. Such flows can add to column densities enough to provide a wind 
optical depth, r w < 1. 

The column density in Model 1 increases from N co \ = 9.6 ■ 10 18 cm” 2 , at the inclination, 
6 = 45° to IVcoi = 1.49 • 10 20 cm” 2 at 6 = 65°, and N co \ = 6.8 • 10 23 cm” 2 at 6 = 90°. 
The column density in Model 3 increases from N co \ = 8 • 10 18 cm -2 at 0 — 45° to iV co i = 
1.3 ■ 10 2 ° cm” 2 at 9 = 65°, and lV col = 7.4 • 10 23 cm” 2 at 6 = 90°. The column densities 
for the denser Models 5-9 (tt = 1.8) are higher. For example, in Model 5 one obtains: 
lV co i = 8.8 • 10 18 cm” 2 at 6 = 45°; N co \ = 3 • 10 20 cm” 2 at 0 = 65°; and N co \ = 2.3 • 10 24 cm” 2 
at 6 = 90°. 

As the luminosity is reduced by almost an order of magnitude, such as in Model 8, part 
of the wind closer to the disk (higher inclinations) becomes denser, and the column density 
increases correspondingly, to N co \ = 8.2 • 10 21 cm” 2 at 0 — 45°; N co i = 3.8 • 10 22 cm” 2 at 
0 = 65°; and N col = 2.4 • 10 24 cm' 2 at 0 = 90°. 

In order to provide obscuration, for example, at 30° away from a disk plane the AGN 
torus should have equatorial densities at least of the order of 10 6 cm” 3 

We also calculated the infrared (i.e. with respect to dust) optical depth, ttr^ in z- 
direction as measured from the upper boundary of the domain. Models 5-9 undergo the 
transition from optically thick to marginally optically thin ones at z ~ 0.25 — 0.5 with a 
shape which closely follows density distribution. The inner low density funnel is always 
optically thin. Above this region there is an extended marginally optically thick region. 
Lower density models 1-4 have optically thick but geometrically thin disk, of of vertical 
extent 0 < 0.2. The rest of the wind has tir^ < 1 closely tracing the distribution of p. 

The distribution of v R (e.g. Figure [2]), demonstrates that v R is increasing in a quasi- 
monotonic way along the spherical radius, r. This is analogous to a one- dimensional, 
radiation-driven stellar wind. On the other hand, the acceleration region is clearly seen 
in the plots of v z . At a given height, z the region of a rapid increase of v z extends in radial 
direction until approximately the optical depth in the infrared (as measured from the right 
boundary), tir,# > 1. When Tjr^ drops, the radiation flux becomes free streaming and no 
significant lifting force in z -direction is generated. 
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5. 0. 5. The mass of the torus 


As the characteristic density, no increases by a more than an order of magnitude when 
we go from Models 1-4 to Models 5-9, the mass of the torus, M tor increases only by a factor 
of 3: from 1.3 x 1O 4 M 0 (Models 1-4), to 4.5 x 1O 4 M 0 (Models 5-9). These numbers provide 
probably the first self-consistent estimates of the torus mass. From the above one can see 
that, though M tor -C Mbh, the self- gravity may be important inside the torus body closer 
to the equatorial plane where densities are higher. We performed test calculations and 
found that basic parameters of the torus are almost unaffected by the choice of the density 


distribution parameter, d in (30). The more important driving parameter which sets the 


scaling for mass-loss rate, velocity field, and the torus mass is the characteristic density no 
(or t t ). 


5.0.6. Relevance to previous work and limitations of the present model 

It is instructive to contrast the approach of the current work to that of Paper I. Com- 
pared to the current work, Paper I contained several simplifying assumptions, from which 
the most stringent were the following two: a) A monotonically accelerating wind with v~ > 0 
was assumed to exist everywhere in the domain and consequently no inflow solutions were 
permitted, and b) only the 0 -component of the velocity was taken into account. In Paper 
I the 2D distribution of E was calculated from the condition V • F = 0 assuming t T <C tf 
(see the discussion in Section [2]). Additionally, in Paper I the equations included only the 
momentum and continuity equation for the matter, and no gas pressure nor the energy equa- 
tion for the matter were considered. Essentially, in the previous studies the radiation force 
was calculated from a solution of a 2D diffusion problem assuming radiative equilibrium, 
and only after that plugging a 0 -component of this force into a ID wind problem in the 
^-direction. 

In this paper we solve the full system of radiation-hydrodynamics equations, adopting 
from Paper I only the description of the boundary conditions. In the present work v max 
is systematically larger than v mSLX = u 0max from Paper I. In Paper I the streamlines were 
directed along the 0 -axis, but in 2D they can bend, return or be tightly packed towards the 
equatorial plane. Thus the addition of another degree of freedom in the present work results 
in smaller mass-loss rates. 

For example, comparing models with T = 0.8, r T ~ 2, we see that w max is considerably 
larger in full RHD modeling: ~ 160kms _1 versus u max — 760kms _1 where results 

from Paper I are marked by a ”*”. This illustrates our previous statements about the very 
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low density and high velocity portion of the wind. The velocities in the denser parts are in 
accord with previous results: (v) ~ 112 kms -1 . Removing the restrictions of Paper I, the 
most stringent being the z-only motion of the gas, also results in higher kinetic energy of 
the wind: Akin — 5 x 10 41 ergs -1 versus L£ in ~ 5.5 x 10 39 ergs -1 . In both sets of models Akin 
is still a tiny fraction of Tboi- The mass-loss rate in the full case is lower than in Paper I: 
for the same set of parameters as above M is reduced from 9.5 M 0 yr -1 to 5 M 0 yr -1 . Lower 
density models have a similar trend of increasing the u max and A y j n compared to the previous 
studies and reducing M . For the models with T = 0.3, tt — 0.5 the latter is again reduced 
by almost 50%, from 2.76 M 0 yr -1 to 1.38 M 0 yr -1 . 


One of the serious limitations of our model is due to the non-vanishing dust opacity. In 
more realistic simulations the very low density parts of the wind, i.e. the funnel seen in the 
p plots will probably not exist. The dust there will most likely sublime and the funnel will 
be filled with much hotter gas. One can expect a picture will resemble an X-ray evaporative 


flow of Dorodnitsyn et al. (2008) but instead of a quasi-static torus there will be an X- 


ray induced evaporation of a dense infrared driven flow. Even without X-rays in the very 
low density part the dust will decouple from the gas, breaking the one-fluid hydrodynamics 
approximation, and thus such dust even if survived will be quickly blown away. Summing 
up we believe that the high velocity and low density component is most likely an artifact of 
our simplifying assumptions. We will address this in future work. 

The finite optical depth provided by some other gas/wind at smaller radii is also not 
taken into account. For example, if a broad absorption line (BAL) wind is formed closer to 
BH and the transverse (approximately) optical depth, r tr between the corona and the dusty, 
infrared-dominated outflow is large, then nothing will be left for the torus. On the other 
hand to have r tr > 1. requires a wind which is more massive than we typically observe in 
BALs. In any case if r tr 1. the model developed in this paper is not applicable. In case 
of no external heating, the only source of radiation is the accretion disk itself. This should 
be addressed in a future work. Help may come from hard X-rays which can penetrate much 


deeper in the torus body providing distributed sources of heat. That was shown by Chang 


et al. (2007) to be quite effective in puffing up the initially geometrically thin accretion 


disk. In more complete simulations, when heating by hard X-ray is implemented we expect 
more efficient acceleration at larger R. At present our calculation gives the most conservative 
estimate of the radiation driving at large distances from the BH. The problem of AGN winds 
is very non-linear and interdependent one, and to answer these questions global multi-group 
radiation-hydrodynamics simulations of AGN accretion disk + winds are required. 
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6. Observational consequences and conclusions 


It is widely accepted that AGN unification schemes require an obscuring toroidal struc- 
ture as a basic premise. Obscuration can be quasi-static or dynamical (winds). An example 
of the first kind of obscuration is the polytropic torus, which is a rotationally supported 
torus with the equation of state P g ~ p 1+1 / ri , where n is the index of the polytrope. Such a 
torus has been shown to be unstable to 3D non-axisymmetric perturbations by |Papaloizou 


& Pringle (1984), but even without such a difficulty it is likely inappropriate as an AGN 


torus prototype. The temperature in a gas-pressure-supported torus is of the order of the 
virial temperature which is 2.6 x 10 6 M 7 /r pc K. That is far too high to be reconciled with 
the existence of dust which requires temperatures of the order of 100 — 1000 K. Obscuration 
can be clumpy/cloudy but we believe that this is the level of complexity which is of the next 
order compared to the fundamental question approached in the current paper. 


The problem of AGN unification via toroidal obscuration can be formulated in the 
following form: what supports the torus against vertical collapse to a geometrically thin 
state and thus maintains its aspect ratio h/R ~ 1? Infrared pressure on dust grains seems 
to be the best candidate. In Paper I it was shown that if the temperature inside the torus 
is of the order of T v ir,r — 312 (n 5 M 7 /rp C ) 1//4 — 987 {n 7 M 7 /r vc ) l ^K an equilibrium between 
radiation pressure, rotational and gravitational forces cannot be maintained, which results 
in formation of an outflow resembling that of an accretion disk wind. 


The conversion of external UV and X-ray radiation into IR pumps the torus with IR 
photons. Internal temperatures of a fewx 10 2 K provide conditions for the extensive presence 
of dust which results in a strong coupling between the gas and IR photon field due to 
the high opacity of dust to IR radiation. In the bulk of the infrared supported torus the 
infrared radiation pressure II P g , with the possible exception of the equatorial region. 
Correspondingly, the torus problem must be formulated in terms of equations of radiation 
hydrodynamics. In the current paper we solved a full system of such equations and found 
that if external BH luminosity exceeds ~ 0.05 L e dd an outflow is created, and that significant 
masses are involved in global motions. 

The obtained mass- loss rates depend on V = L/L edd and on the characteristic density, 
no which scales the distribution of p in the equatorial plane. Generally, models with higher V 
and with larger no (or alternatively tt) tend to have higher mass-loss rates. However there 
exist an overlap, when marginally optically thin but luminous models produce mass-loss 
rates similar to the optically thick and less luminous ones. For example, the model with 
tt — 0.5, and T = 0.3 has mass-loss rate, M = 1.4M 0 yr _1 , similar to the model with 
tt — 2 and T = 0.05. Besides such overlap, in optically-thick models with higher T much 
larger M is obtained than in optically thin ones. For example, in one of the Thomson thick 
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and luminous models with tt — 2 and F = 0.8 the mass-loss rate is M = 5 M 0 yr^ 1 Another 
important finding is that most of the flow has velocity (v) < f/ esc , that is too small to escape 
the potential well of the BH. Thus the AGN obscuring flow is better described in terms of a 
failed wind rather than viewed as a bipolar and 2.5D analogy of a stellar wind. 


In our calculations a number of simplifications were made. The most restrictive one is 
the shape of the conversion layer, where UV and soft X-rays are converted into the infrared. 
We assume that such a conversion is happening at the vertical boundary of the computational 
domain while in real AGN two important complications arise: i) the curvature of such layer is 
important and can only be determined during self-consistent global radiation hydrodynamics 
simulations, and ii) the corresponding ’’photospheres” are located at different depths which 
could be captured in RHD simulations with multi-group description of radiation. Taking 
into account additional heating by hard X-rays would also contribute to the structure and 
dynamics of the obscuring flow (Chang et al. 2007 Shi & Krolik 2008). These effects must 
be incorporated into future global radiation hydrodynamics simulations. 


Observationally, tracing wind kinematics through the detection of maser emission is one 
way to ’’see” an AGN obscuring flow. From our results it follows that infrared supported flow 
naturally produces outflows with bulk velocities as large as ~ few x lOOkms -1 . It is beyond 
the scope of this paper to calculate conditions and particular locations in the flow suitable for 
maser emission. However, our solutions allow us to predict that if such emission is observed 
at distances 0.4 — 1.5 from a BH and being offset from the corresponding Keplerian velocity 
of the equatorial disk by several hundreds kms" 1 , it may indicate an IR-driven outflow. 
Such evidence may already be present in the VLBI observations of a broad bipolar outflow 
in Circinus galaxy in H 2 0 maser emission (Greenhill et al. 2003). 
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